"""
Use the Ratio of Uniforms algorithm to sample a distribution
============================================================
"""

# %%
# In this example, we illustrate in dimension :math:`\inputDim = 1` the ratio of uniforms algorithm  to sample a distribution which
# probability density function is known up to a normalization factor.
#
# Consider a scalar distribution whose
# probability density function :math:`p` is known up to a normalization factor:
# let :math:`f: \Rset \rightarrow \Rset^+` be a function such that
# :math:`p(x) = cf(x)` for all :math:`x \in \Rset` with :math:`c \in \Rset^+_*`.
#
# Let :math:`r \in \Rset^+_*`. Let  :math:`A_{f,r}` be the domain defined by:
#
# .. math::
#
#    A_{f,r} = \left\{ (u, v) \in \Rset^2 \, |\, 0 \leq u \leq f\left(\dfrac{v}{u^r}\right)^{\frac{1}{1+r}}\right \}
#
# The Lebesgue measure of :math:`A_{f,r}`  is equal to :math:`\dfrac{1}{c(1+r)}`.
#
# Let :math:`(U, V)` be a random variable uniformly distributed on :math:`A_{f,r}`.
# Then, :math:`X = \dfrac{V}{U^r}` is a random variable on
# :math:`\Rset` distributed according to :math:`p`.
#
# Under the condition that :math:`f(x)^{\frac{1}{1+r}}` and
# :math:`x f(x)^{\frac{1}{1+r}}` are bounded, then the domain
# :math:`A_{f,r}` is bounded and can be included in a rectangular bounding box :math:`\tilde{A}_{f,r}`
# defined by:
#
# .. math::
#
#    \tilde{A}_{f,r} & = \left[0, \sup_{x} f(x)^{\frac{1}{1+r}} \right] \times
#    \left[ \inf_{x} x f(x)^{\frac{1}{1+r}} ,  \sup_{x} x f(x)^{\frac{1}{1+r}}\right]\\
#                    & = \left[0, u_{sup} \right] \times \left[v_{inf}, v_{sup}\right]
#
# Note that :math:`u_{sup} > 0`, :math:`v_{sup} \geq 0` and  :math:`v_{inf} \leq 0`.
#
# This allows one to sample uniformly the domain :math:`A_{f,r}` by rejection sampling inside
# :math:`\tilde{A}_{f,r}`.
#
# Study of the :math:`(u,v)`-space.
# ---------------------------------
# The boundary of :math:`\tilde{A}_{f,r}` is defined by:
#
# .. math::
#
#     \partial \tilde{A}_{f,r} = \cF_1 \cup \cF_2 \cup \cF_3 \cup \cF_4
#
# where the boundaries :math:`\cF_i` are defined by:
#
# .. math::
#
#    \cF_1 & = \left \{ (u, v_{sup}), 0 \leq u \leq u_{sup} \right \} \\
#    \cF_2 & = \left \{ (u, v_{inf}), 0 \leq u \leq u_{sup} \right \}\\
#    \cF_3 & = \left \{ (u_{sup}, v), v_{inf} \leq v \leq v_{sup} \right \}\\
#    \cF_4 & = \left \{ (0, v), v_{inf} \leq v \leq v_{sup} \right \}
#
# The boundary of :math:`A_{f,r}` is defined by:
#
# .. math::
#
#    \partial A_{f,r}  = \left\{ (u,v) = \left( f(x)^{\dfrac{1}{1+r}}, x f(x)^{\dfrac{r}{1+r}} \right),
#    x \in \Rset \right\}
#
# The following function draws the domain :math:`A_{f,r}` and the bounding box :math:`\partial \tilde{A}_{f,r}`
# in the :math:`(u,v)`-space.
import openturns as ot
import openturns.viewer as otv
import openturns.experimental as otexp
from math import pi, exp


def draw_RatioOfUniforms_UV(ratio_algo):
    r = ratio_algo.getR()
    range_logF = ratio_algo.getRange()
    xInf = range_logF.getLowerBound()[0]
    xSup = range_logF.getUpperBound()[0]
    uSup = ratio_algo.getSupU()
    vInf = ratio_algo.getInfV()[0]
    vSup = ratio_algo.getSupV()[0]
    logF = ratio_algo.getLogUnscaledPDF()

    # \partial A_{f,r} domain
    x = xInf
    u = 0.0
    v = 0.0
    xPrec = None
    uPrec = None
    vPrec = None
    first = True
    step = 0.01 * (uSup + vSup - vInf)
    delta = step
    data_UV = ot.Sample(0, 2)
    while x < xSup:
        val_F = exp(logF([x])[0])
        u = val_F ** (1.0 / (1.0 + r))
        v = x * val_F ** (r / (1.0 + r))
        # Step adaptation after the first point
        if not first:
            # If the step is too large
            while abs(u - uPrec) + abs(v - vPrec) > 2 * delta:
                step *= 0.5
                x = xPrec + step
                if x < xSup:
                    val_F = exp(logF([x])[0])
                    u = val_F ** (1.0 / (1.0 + r))
                    v = x * val_F ** (r / (1.0 + r))
            # If the step is too small
            while abs(u - uPrec) + abs(v - vPrec) < 0.5 * delta and x < xSup:
                step *= 2.0
                x = xPrec + step
                if x < xSup:
                    val_F = exp(logF([x])[0])
                    u = val_F ** (1.0 / (1.0 + r))
                    v = x * val_F ** (r / (1.0 + r))
        xPrec = x
        uPrec = u
        vPrec = v
        x += step
        first = False
        data_UV.add([u, v])
    # to get a closed boundary
    data_UV.add(data_UV[0])
    g = ot.Graph("", "u", "v", True, "bottomright")
    data_A = ot.Curve(data_UV)
    data_A.setColor("purple")
    g.add(data_A)

    # \partial \tilde{A}_{f,r} bounding box
    boundary_1 = ot.Curve([[0, vSup], [uSup, vSup]])
    boundary_1.setColor("blue")
    boundary_2 = ot.Curve([[0, vInf], [uSup, vInf]])
    boundary_2.setColor("red")
    boundary_3 = ot.Curve([[uSup, vInf], [uSup, vSup]])
    boundary_3.setColor("green")
    boundary_4 = ot.Curve([[0, vInf], [0, vSup]])
    boundary_4.setColor("black")
    g.add(boundary_1)
    g.add(boundary_2)
    g.add(boundary_3)
    g.add(boundary_4)

    leg = g.getLegends()
    leg[0] = r"boundary of $A_{f,r}$"
    leg[1] = r"$\mathcal{F}_1$"
    leg[2] = r"$\mathcal{F}_2$"
    leg[3] = r"$\mathcal{F}_3$"
    leg[4] = r"$\mathcal{F}_4$"
    g.setLegends(leg)

    return g


# %%
# Study of the :math:`(x,y)`-space.
# ---------------------------------
# We use the mapping function :math:`\varphi` defined by:
#
# .. math::
#
#    \varphi: & \Rset^2 \rightarrow \Rset^2 \\
#             & (u,v)  \rightarrow \left( \dfrac{v}{u^r}, u \right)
#
# Then we have:
#
# .. math::
#
#    \varphi(A_{f,r}) = \left\{ (x,y) \, |\, x \in \Rset, 0 \leq y \leq f(x)^{\dfrac{1}{1+r}} \right\}
#
# and:
#
# .. math::
#
#     \partial \varphi(\tilde{A}_{f,r}) = \varphi(\cF_1) \cup \varphi(\cF_2) \cup \varphi(\cF_3)
#
# where the boundaries :math:`\varphi(\cF_i)` are defined by:
#
# .. math::
#
#    \varphi(\cF_1) & = \left\{ \varphi(u, v_{sup}), 0 \leq u \leq u_{sup} \right\}\\
#          & = \left\{ \left(x,\left(\dfrac{v_{sup}}{x}\right)^{1/r} \right),
#              x \geq \dfrac{v_{sup}}{u_{sup}^{r}} \right \}
#
# .. math::
#
#    \varphi(\cF_2) & = \left \{ \varphi(u, v_{inf}), 0 \leq u \leq u_{sup} \right\}\\
#          & = \left\{ \left(x,\left(\dfrac{v_{inf}}{x}\right)^{1/r}\right),
#              x \leq \dfrac{v_{inf}}{u_{sup}^{r}} \right \}
#
# .. math::
#
#    \varphi(\cF_3) & = \left\{ \varphi(u_{sup}, v), v_{inf} \leq v \leq v_{sup} \right\}\\
#          & = \left\{ (x, u_{sup}),
#              \dfrac{v_{inf}}{u_{sup}^{r}} \leq x \leq \dfrac{v_{sup}}{u_{sup}^{r}}  \right\}
#
# The following function draws the boundary :math:`\partial \varphi(A_{f,r})` and the bounding box
# :math:`\partial \tilde{A}_{f,r}` in the :math:`(x,y)`-space.


def draw_RatioOfUniforms_XY(ratio_algo, nPoints):
    r = ratio_algo.getR()
    range_logF = ratio_algo.getRange()
    xInf = range_logF.getLowerBound()[0]
    xSup = range_logF.getUpperBound()[0]
    uSup = ratio_algo.getSupU()
    vInf = ratio_algo.getInfV()[0]
    vSup = ratio_algo.getSupV()[0]
    logF = ratio_algo.getLogUnscaledPDF()

    # \partial varpĥi(A_{f,r}) domain
    eps = 1e-10
    step = (xSup - xInf - 2.0 * eps) / (nPoints - 1)
    data_XY = ot.Sample(nPoints, 2)
    for i in range(nPoints):
        x = xInf + eps + i * step
        val_F = exp(logF([x])[0])
        data_XY[i, 0] = x
        data_XY[i, 1] = val_F ** (1.0 / (1.0 + r))
    g = ot.Graph("", "x", "y", True, "lower center")
    data_A = ot.Curve(data_XY)
    data_A.setColor("purple")
    g.add(data_A)

    # \partial varphi(\tilde{A}_{f,r}) bounding box
    leg = ot.Description()
    leg.add(r"boundary of $\varphi(A_{f,r})$")
    boundary_1_fct = ot.SymbolicFunction(
        "x", "(" + str(vSup) + "/x)^(1/" + str(r) + ")"
    )
    if vSup / (uSup**r) < xSup:
        boundary_1 = boundary_1_fct.draw(vSup / (uSup**r), xSup, 101)
        boundary_1.setColors(["blue"])
        leg.add(r"$\varphi(\mathcal{F}_1)$")
        g.add(boundary_1)
    if vInf / (uSup**r) < uSup:
        boundary_2_fct = ot.SymbolicFunction(
            "x", "(" + str(vInf) + "/x)^(1/" + str(r) + ")"
        )
        boundary_2 = boundary_2_fct.draw(xInf, vInf / (uSup**r), 101)
        boundary_2.setColors(["red"])
        leg.add(r"$\varphi(\mathcal{F}_2)$")
        g.add(boundary_2)
    boundary_3 = ot.Curve([[vInf / (uSup**r), uSup], [vSup / (uSup**r), uSup]])
    boundary_3.setColor("green")
    g.add(boundary_3)
    leg.add(r"$\varphi(\mathcal{F}_3)$")
    g.setLegends(leg)

    return g


# %%
# Illustration.
# -------------
# We create the function :math:`\log f` defined by:
#
# .. math::
#    :label: robertEx
#
#    \log f(x) = \log(\cos(x)) + x,  \quad  x \in \left[-\dfrac{\pi}{2}, \dfrac{\pi}{2}\right]
#
# and the :class:`~openturns.experimental.RatioOfUniforms`.


log_UnscaledPDF = ot.SymbolicFunction("x", "log(cos(x)) + x")
eps = 1e-5
range_PDF = ot.Interval(-pi / 2.0 + eps, pi / 2.0 - eps)
ratio_Algo = otexp.RatioOfUniforms(log_UnscaledPDF, range_PDF)

# %%
# We draw the boundary :math:`\partial A_{f,r}` and the bounding box :math:`\partial \tilde{A}_{f,r}`
# in the :math:`(u,v)`-space.
g_UV = draw_RatioOfUniforms_UV(ratio_Algo)
g_UV.setTitle(
    rf"$r = ${ratio_Algo.getR():.3g}, $\tau = ${ratio_Algo.getAcceptanceRatio():.3g}"
)

# %%
# We draw the boundary :math:`\partial \varphi(A_{f,r})` and the bounding box :math:`\partial \varphi(\tilde{A}_{f,r})`
# in the :math:`(x,y)`-space.
# sphinx_gallery_thumbnail_number = 1
g_XY = draw_RatioOfUniforms_XY(ratio_Algo, 1001)
g_XY.setTitle(
    rf"$r = ${ratio_Algo.getR():.3g}, $\tau = ${ratio_Algo.getAcceptanceRatio():.3g}"
)

grid_default = ot.GridLayout(1, 2)
grid_default.setGraph(0, 0, g_UV)
grid_default.setGraph(0, 1, g_XY)
grid_default.setTitle(r"$\log f(x) = \log(\cos(x)) + x, x \in ]-\pi/2, \pi/2[$")
view = otv.View(grid_default)

# %%
# We change the :math:`r` value and we draw the same plots.
grid_full = ot.GridLayout(4, 2)
for i, r in enumerate([0.5, 0.8, 1.2, 2.0]):
    ratio_Algo.setR(r)
    g_UV = draw_RatioOfUniforms_UV(ratio_Algo)
    g_UV.setLegendPosition("")
    g_UV.setTitle(
        rf"$r = ${ratio_Algo.getR():.3g}, $\tau = ${ratio_Algo.getAcceptanceRatio():.3g}"
    )

    g_XY = draw_RatioOfUniforms_XY(ratio_Algo, 1001)
    g_XY.setLegendPosition("")
    g_XY.setTitle(
        rf"$r = ${ratio_Algo.getR():.3g}, $\tau = ${ratio_Algo.getAcceptanceRatio():.3g}"
    )

    grid_full.setGraph(i, 0, g_UV)
    grid_full.setGraph(i, 1, g_XY)

grid_full.setTitle(r"$\log f(x) = \log(\cos(x)) + x,  \quad x \in ]-\pi/2, \pi/2[$")
view = otv.View(grid_full)

# %%
# Study of the acceptance ratio.
# ------------------------------
# We study the acceptance ratio :math:`\tau(r)` as a function of the parameter :math:`r` for some distributions.
# In these cases, the value is exact and computed from the volume of :math:`A_{f,r}` and the one of :math:`\tilde{A}_{f,r}`.
# We draw the function :math:`r \rightarrow \tau(r)`. We conclude that the default value :math:`r=1` seems to be correct
# in all these cases.
#
mixt_gaus = ot.Mixture([ot.Normal(-2.0, 1.0), ot.Normal(4.0, 1.0)], [0.2, 0.8])
logPDF_list = [
    ot.Normal().getLogPDF(),
    ot.Normal(2).getLogPDF(),
    ot.Normal(3).getLogPDF(),
    ot.Exponential().getLogPDF(),
    mixt_gaus.getLogPDF(),
]
range_list = [
    ot.Normal().getRange(),
    ot.Normal(2).getRange(),
    ot.Normal(3).getRange(),
    ot.Exponential().getRange(),
    mixt_gaus.getRange(),
]
name_list = ["Normal(1)", "Normal(2)", "Normal(3)", "Exponential()", "Mixture Normal"]

g = ot.Graph(
    r"Acceptance ratio $\tau$ as a function of $r$", r"$r$", r"$\tau$", True, "topright"
)
g.setLogScale(1)
numberR = 50
a = (1e2) ** (2.0 / numberR)
data_r_ratio = ot.Sample(numberR + 1, 2)
isScaled = True
for k in range(len(logPDF_list)):
    ratio_Algo = otexp.RatioOfUniforms(logPDF_list[k], range_list[k], isScaled)
    for i in range(numberR + 1):
        r = a ** (i - 0.5 * numberR)
        ratio_Algo.setR(r)
        accept_Ratio = ratio_Algo.getAcceptanceRatio()
        data_r_ratio[i] = [r, accept_Ratio]
    draw = ot.Curve(data_r_ratio)
    draw.setLegend(name_list[k])
    g.add(draw)

isScaled = False
ratio_Algo = otexp.RatioOfUniforms(log_UnscaledPDF, range_PDF, isScaled)
for i in range(numberR + 1):
    r = a ** (i - 0.5 * numberR)
    ratio_Algo.setR(r)
    accept_Ratio = ratio_Algo.getAcceptanceRatio()
    data_r_ratio[i] = [r, accept_Ratio]
draw = ot.Curve(data_r_ratio)
draw.setLegend(
    r"$f(x) = \cos(x)e^x\mathbf{1}_{\left[-\frac{\pi}{2}, \frac{\pi}{2}\right]}(x)$"
)
g.add(draw)

g.setLegendPosition("upper left")
g.setLegendCorner([1.0, 1.0])
view = otv.View(g)

# %%
# Exploration of the optimization problems.
# -----------------------------------------
# We illustrate how to explore the results of the optimization problems that have been solved to define
# :math:`\tilde{A}_{f,r}`.
# We consider the example defined by the function:
#
# .. math::
#
#     f(x) = \frac{1}{2} (2 + \sin(x)^2) \exp \left[- \left(2 + \cos(3x)^3 + \sin(2x)^3 \right) x
#     \right]  \mathbf{1}_{[0, 2 \pi]}(x).
#
# The method *initialize* returns all the starting points that have been used to compute respectively :math:`u_{sup}`,
# :math:`v_{inf}` and :math:`v_{sup}`.
# The maximum number of starting points is defined in :class:`~openturns.ResourceMap`, entry
# *RatioOfUniformas-MaximumMultiStart*. These
# points are selected scrolling through  the Sobol low discrepance sequence in dimension :math:`\inputDim =  1`
# and which maximum size is defined in :class:`~openturns.ResourceMap`, entry *RatioOfUniforms-CandidateNumber*.
f = ot.SymbolicFunction(
    "x", "0.5 * (2 + sin(x)^2) * exp( -( 2 + cos(3*x)^3 + sin(2*x)^3) * x )"
)
log_UnscaledPDF = ot.ComposedFunction(ot.SymbolicFunction("x", "log(x)"), f)
lower_bound = 0.0
upper_bound = 2.0 * pi
range_PDF = ot.Interval(lower_bound, upper_bound)

graph = f.draw(lower_bound, upper_bound, 512)
graph.setTitle(
    r"Christian Robert function: $f(x) =  0.5(2 + \sin^2 x) e^{ -x( 2 + \cos^3 3x + \sin^3 2x)}, x \in [0, 2\pi]$"
)
graph.setXTitle(r"$x$")
graph.setYTitle(r"$f(x)$")
view = otv.View(graph)

ratio_Algo = otexp.RatioOfUniforms(log_UnscaledPDF, range_PDF)

# %%
# In order to get all the intermediate results, we change the :class:`~openturns.ResourceMap`,
# entry *MultiStart-KeepResults*.
#
# Furthermore, as :math:`I = [0, 2 \pi]=  [a_i, b_i]`, then :math:`a_i \geq 0` implies that :math:`v_{inf} = 0`.
# There are only the bounds :math:`u_{sup}` and :math:`v_{sup}` to find: the collection of optimization problems
# only contains 2 items.
ot.ResourceMap.SetAsBool("MultiStart-KeepResults", True)
coll_Multi_Start = ratio_Algo.initialize()
start_Sample_SupU = coll_Multi_Start[0].getStartingSample()
start_Sample_SupV = coll_Multi_Start[1].getStartingSample()

# %%
# We plot on the previous graph all the starting points (red) used to find :math:`u_{sup}` and the associated maxima found by the algorithm. We
# also link each starting point to its associated optimal value (green).
Cloud_start_Usup = ot.Cloud(
    start_Sample_SupU, ot.Sample(start_Sample_SupU.getSize(), [0.0])
)
Cloud_start_Usup.setColor("red")
Cloud_start_Usup.setPointStyle("fsquare")
graph.add(Cloud_start_Usup)

coll_results = coll_Multi_Start[0].getResultCollection()
sample_optima = ot.Sample(0, 2)
for i in range(coll_results.getSize()):
    opt_x = coll_results[i].getOptimalPoint()[0]
    opt_value = f([opt_x])[0]
    sample_optima.add([opt_x, opt_value])
    line_convergence = ot.Curve([[start_Sample_SupU[i, 0], 0.0], [opt_x, opt_value]])
    line_convergence.setLineStyle("dashed")
    line_convergence.setColor("black")
    graph.add(line_convergence)
cloud_optima = ot.Cloud(sample_optima)
cloud_optima.setColor("green")
cloud_optima.setPointStyle("star")
graph.add(cloud_optima)
graph.setYTitle("")
view = otv.View(graph)


# %%
view.ShowAll()
